﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace SmartMathLibrary.FFTPack
{
    /// <summary>
    /// FFT transform of a real periodic sequence. 
    /// </summary>
    [Serializable]
    public class RealDoubleFFT : RealDoubleFFTMixed
    {
        /// <summary>
        /// This field holds the norm factor.
        /// </summary>
        private double normFactor;

        /// <summary>
        /// This field holds the wave table.
        /// </summary>
        private readonly double[] wavetable;

        /// <summary>
        /// This field holds the dimension.
        /// </summary>
        private readonly int ndim;

        /// <summary>
        /// Initializes a new instance of the <see cref="RealDoubleFFT"/> class.
        /// </summary>
        /// <param name="n">The size of a real data sequence. When n is a multiplication of small 
        /// numbers (4, 2, 3, 5), this FFT transform is very efficient.</param>
        public RealDoubleFFT(int n)
        {
            if (15 + 2 * n < 0)
            {
                throw new ArgumentException("15 + 2 * n < 0", "n");
            }

            ndim = n;
            normFactor = n;

            if (wavetable == null || wavetable.Length != (2 * ndim + 15))
            {
                wavetable = new double[2 * ndim + 15];
            }

            Rffti(ndim, wavetable);
        }

        /// <summary>
        /// This field can be used to normalize this FFT transform. This is because a call of forward 
        /// transform followed by a call of backward transform will multiply the input sequence by the 
        /// norm factor.
        /// </summary>
        /// <value>The norm factor to use for the transform.</value>
        public double NormFactor
        {
            get { return normFactor; }
            set { normFactor = value; }
        }

        /// <summary>
        /// Forward real FFT transform. It computes the discrete transform of a real data sequence.
        /// </summary>
        /// <param name="x">An array which contains the sequence to be transformed. After FFT, x contains the 
        /// transform coeffients used to construct n complex FFT coeffients.</param>
        public void ForwardTransform(double[] x)
        {
            if (x == (double[]) null)
            {
                throw new ArgumentNullException("x");
            }

            if (x.Length != ndim)
            {
                throw new ArgumentException("The length of data can not match that of the wavetable");
            }

            Rfftf(ndim, x, wavetable);
        }

        /// <summary>
        /// Forward real FFT transform. It computes the discrete transform of a real data sequence.
        /// </summary>
        /// <param name="x">An array which contains the sequence to be transformed. After FFT, x contains the 
        /// transform coeffients used to construct n complex FFT coeffients.</param>
        /// <param name="y">The first complex n/2 (when n is odd) or (n/2+1) (when n is even) FFT coeffients.</param>
        public void ForwardTransform(double[] x, Complex1D y)
        {
            if (x.Length != ndim)
            {
                throw new ArgumentException("The length of data can not match that of the wavetable");
            }
            Rfftf(ndim, x, wavetable);

            if (ndim % 2 == 0)
            {
                y.X = new double[ndim / 2 + 1];
                y.Y = new double[ndim / 2 + 1];
            }
            else
            {
                y.X = new double[(ndim + 1) / 2];
                y.Y = new double[(ndim + 1) / 2];
            }


            y.X[0] = x[0];
            y.Y[0] = 0.0D;
            for (int i = 1; i < (ndim + 1) / 2; i++)
            {
                y.X[i] = x[2 * i - 1];
                y.Y[i] = x[2 * i];
            }
            if (ndim % 2 == 0)
            {
                y.X[ndim / 2] = x[ndim - 1];
                y.Y[ndim / 2] = 0.0D;
            }
        }

        /// <summary>
        /// Backward real FFT transform. It is the unnormalized inverse transform of ForwardTransform(double[]).
        /// </summary>
        /// <param name="x">An array which contains the sequence to be transformed. After FFT, x contains 
        /// the transform coeffients. Also see the comments of ForwardTrasform()(double[]).</param>
        public void BackTransform(double[] x)
        {
            if (x == (double[]) null)
            {
                throw new ArgumentNullException("x");
            }

            if (x.Length != ndim)
            {
                throw new ArgumentException("The length of data can not match that of the wavetable");
            }
            Rfftb(ndim, x, wavetable);
        }

        /// <summary>
        /// Backward real FFT transform. It is the unnormalized inverse transform of ForwardTransform(Complex1D, double[]).
        /// </summary>
        /// <param name="x">An array which contains the sequence to be transformed. When n is odd, 
        /// it contains the first n/2 complex data; when n is even, it contains (n/2+1) complex data.</param>
        /// <param name="y">The real FFT coeffients.</param>
        public void BackTransform(Complex1D x, double[] y)
        {
            if (x.X == (double[]) null)
            {
                throw new ArgumentNullException("x");
            }

            if (x.Y == (double[]) null)
            {
                throw new ArgumentNullException("x");
            }

            if (ndim % 2 == 0)
            {
                if (x.X.Length != ndim / 2 + 1)
                {
                    throw new ArgumentException("The length of data can not match that of the wavetable");
                }
            }
            else
            {
                if (x.X.Length != (ndim + 1) / 2)
                {
                    throw new ArgumentException("The length of data can not match that of the wavetable");
                }
            }

            y[0] = x.X[0];
            for (int i = 1; i < (ndim + 1) / 2; i++)
            {
                y[2 * i - 1] = x.X[i];
                y[2 * i] = x.Y[i];
            }
            if (ndim % 2 == 0)
            {
                y[ndim - 1] = x.X[ndim / 2];
            }
            Rfftb(ndim, y, wavetable);
        }
    }
}